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We extend previous work on high energy nuclear collisions in the Color Glass Condensate model 
to study collisions of finite ultrarelativistic nuclei. The changes implemented include a) imposition 
of color neutrality at the nucleon level and b) realistic nuclear matter distributions of finite nuclei. 
The saturation scale characterizing the fields of color charge is explicitly position dependent, A s = 
A s (xt). We compute gluon distributions both before and after the collisions. The gluon distribution 
in the nuclear wavefunction before the collision is significantly suppressed below the saturation scale 
when compared to the simple McLerran- Venugopalan model prediction, while the behavior at large 
momentum pr A s remains unchanged. We study the centrality dependence of produced gluons 
and compare it to the centrality dependence of charged hadrons exhibited by the RHIC data. We 
demonstrate the geometrical scaling property of the initial gluon transverse momentum distributions 
for different centralities. Classical Yang-Mills results for pr < A s are simply matched to perturbative 
QCD computations for pr > A 3 -the resulting energy per particle is significantly lower than the purely 
classical estimates. Our results for nuclear collisions can be used as initial conditions for quantitative 
studies of the further evolution and possible equilibration of hot and dense gluonic matter produced 
in heavy ion collisions. Finally, we study pA collisions within the classical framework. Our results 
agree well with previously derived analytical results in the appropriate kinematical regions. 



PACS numbers: 24.85.+p,25.75.-q,12.38.Mh 

I. INTRODUCTION 

An outstanding problem of considerable theoretical 
and experimental interest is the possible formation of a 
quark gluon plasma (QGP) in collisions of nuclei at very 
high energies. Since the collisions are very violent and 
the time scales involved are ephemeral, the answer to the 
question of whether a QGP is formed depends sensitively 
on the initial conditions for the matter produced in the 
collision. Specifically, it requires an understanding of the 
distributions of partons in the wavefunctions of the two 
nuclei before the collision. What do these parton distri- 
butions look like? 

It is well known that multi-particle production in high 
energy collisions is dominated by modes in the nuclear 
wavefunction which carry a small momentum fraction x 
of the nuclear momentum Understanding the cor- 
rect initial conditions therefore requires that we under- 
stand the properties of these small x modes. In re- 
cent years, a sophisticated effective field theory approach 
has been developed to describe properties of partons at 
small x 0, Jl E|-these form a Color Glass Condensate 
(CGC) g § |T The CGC is characterized by a bulk 
scale A s which grows with energy and the size of the nu- 
clei. For RHIC energies, if one assumes A s is a constant 
as for cylindrical nuclei, then A s - 1-2 GeV @, § |, fl0|. 
The field strengths in the saturation region behave as 
~ l/as'. since as(A s ) <C 1, the field strengths are large. 
Furthermore, the occupation number of saturated gluons 
is also ~ l/as 1< 

The initial conditions for nuclear collisions can be for- 



mulated in the CGC 0, ||, 0] . Since the occupation 
numbers of the partons before the collision are large, the 
very initial stage of the nuclear collision can be treated 
classically. In practice, this means that one solves the 
Yang-Mills equations for two sources of color charge with 
initial conditions given by the classical fields of the two 
nuclei before the collision. This program is difficult to ac- 
complish analytically, but it can be formulated and solved 
numerically ^| . An important assumption is 

that of boost invariance: the gauge fields are assumed not 
to depend on the space-time rapidity r/. The QCD Hamil- 
tonian (in A T — gauge, where r is the proper time) can 
then be formulated as a dimensionally reduced 2+1-D 
theory of transverse gauge fields coupled to an adjoint 
scalar field. Numerically, the lattice analog of the di- 
mensionally reduced Hamiltonian is the Kogut-Susskind 
Hamiltonian [^9) coupled to an adjoint scalar field [ p^| . 
The space-time evolution of the classical gauge fields after 
the collision is obtained by solving, for each configuration 
of the color charge sources, Hamilton's equations for the 
canonically conjugate momenta and fields. At late times 
t > 1 / A s the system is well approximated as a collection 
of weakly coupled harmonic oscillators, and the energy 
and number distributions can be computed accordingly 
for an SU(2) gauge theory |l6|, |l7| and an SU(3) gauge 
theory |jg| after averaging over all color configurations. 

In these earlier studies, nuclear collisions, for simplic- 
ity, were idealized as central collisions of infinite, cylin- 
drical nuclei. The color charge squared A 2 was taken to 
be a constant for the uniform cylindrical nuclei. Fur- 
thermore, color neutrality was imposed only in a global 
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sense ||, namely, the color charge distribution over the 
entire nucleus was constrained to be zero. While very 
useful in obtaining first estimates of the space-time evo- 
lution of the produced gluonic matter, these studies did 
not make predictions for realistic nuclear collisions. In 
addition, studies of the distributions in peripheral colli- 
sions, in particular of the azimuthal anisotropy associ- 
ated with elliptic flow, require finite nuclei and realistic 
nuclear matter distributions within each nucleus. In a re- 
cent paper we implemented all of these realistic require- 
ments to compute elliptic flow in the CGC model [ |l0| . 
In this paper we will discuss the effects of imposing color 
neutrality in more detail and compute energy and num- 
ber distributions as a function of centrality. 

Strictly speaking, the results of our numerical simula- 
tions are only valid for very early times (r ~ 1/A S ) when 
the occupation numbers of the fields are large enough for 
the classical field approximation to hold. The subsequent 
evolution of the system can be treated by solving semi- 
classical transport equations. An interesting scenario of 
how thermalization may be achieved in this picture was 
discussed by Baier et al. [|o[ The results presented 
here can be used as initial conditions for detailed simula- 
tions of the equilibration dynamics in the early stage of 
heavy ion collisions |22|, |24|, |2f| . 

The results here can also be compared to the data on 
global observables from RHIC, if one assumes parton- 
hadron duality || . Remarkably, the CGC results 

on the centrality and energy dependence of charged par- 
ticle multiplicities and on the rapidity and pr distribu- 
tions appear to agree with data. However, there are data 
on the transverse energy per particle and on the elliptic 
flow (i>2 and i^Cpt)) which appear to require additional 
assumptions beyond simple parton-hadron duality. We 
will discuss constraints from the data and from our sim- 
ulations on our scenarios of equilibration and hydrody- 
namic flow in heavy ion collisions. 

Finally, we will discuss pA collisions at RHIC energies. 
This topic has recently attracted much theoretical inter- 
est recently (28| because of the likelihood that such (more 
specifically dA) collisions will be performed at RHIC in 
the immediate future. This problem was first discussed 
in the classical framework by Kovchegov and Mueller p9| ] 
and subsequently by Dumitru and McLerran Both 
these authors derived analytical expressions for gluon 
production in pA collisons. In particular, Dumitru and 
McLerran showed that the gluon distribution, while pro- 
portional to 1/pt for A s i < A S 2 < Pt, is proportional to 
\jp\ in the kinematic region A s i < pt < A S 2- Here A s i 
and A S 2 are the color charges squared of the proton and 
nucleus respectively. Reproducing these results is a test 
of our numerical framework and we will show that our 
results indeed agree well with the analytical expressions. 

This paper is organized as follows. In the following sec- 
tion, we will discuss the numerical simulations of nuclear 
collisions with emphasis on improvements over previous 
work. In particular, we will discuss how the condition of 
color neutrality is imposed on the color charge configu- 



rations to ensure that there are no gluon fields and gluon 
production outside the nucleus. In section HI. we will dis- 
cuss the results of our simulations for gluon distributions 
both before and after the collision. We compute the cen- 
trality dependence of the distribution of produced gluons 



and compare these to the RHIC data. In section [V 



wo 



will discuss the results of our simulation of a pA collision. 
In the final section, we conclude with a summary and a 
discussion of further improvements in our approach. 



II. MODEL OF COLLISIONS OF FINITE 
ULTRARELATIVISTIC NUCLEI 

We will discuss in this section the McLerran- 
Venugopalan model (MV) Q for ultrarelativistic nuclei 
with realistic nuclear matter distributions. We will com- 
pute the classical gluon distributions of the nuclei both 
before and after the collision. 

We first discuss parton distributions in a single nucleus 
in the MV approach [| . In the original MV approach, 
the color charge squared per unit area, A^, is taken to be 
much larger than the confinement scale ( ~ ^qcd)- For 
realistic nuclei at finite energies, confining effects cannot 
be ignored. For instance if in a finite nucleus no color 
neutrality condition is imposed, the gluon distribution 
outside the nucleus would be finite even if the distribu- 
tion of color charges were localized. The issue of color 
neutrality in the MV model of a single nucleus was first 
discussed at length by Lam and Mahlon |nj]. In their 
work, they considered large nuclei with uniform nuclear 
matter distributions. In this section, we will construct 
color charge distributions by imposing different color neu- 
trality constraints at the nucleon level. 

Once we have the color charge distributions, we are 
ready to consider the case of nuclear collisions. The prob- 
lem was first formulated for very large nuclei by Kovner, 
McLerran and Weigert pJj . The full numerical solution 
of the problem for uniform cylindrical nuclei was dis- 
cussed by us in a series of papers (l^, [l6| [l7| We 
will discuss in this section numerical solutions for the 
case of finite nuclei. Color neutrality constraints prove 
essential to avoiding particle production outside interac- 
tion region. 



McLerran Venugopalan model for a finite 
ultrarelativistic nucleus 



We will begin with a very brief review of the McLerran- 
Venugopalan model t|. At high energies, multi- 
particle production is dominated by wee partons which 
carry a small fraction x of the longitudinal light cone mo- 
mentum P + of the nucleus. The McLerran- Venugopalan 
model is a model of the distribution of these wee par- 
tons coupled to hard parton sources at large x. In 
the infinite momentum frame, the light cone current of 
these sources can be expressed as (for a nucleus mov- 
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ing along the positive z-axis) J = (J+, 0,0,0), where 
J + - a = p a (xT)5{x~) is the color charge density on the 
light cone. Note that light cone co-ordinates are defined 
here as x ± = (t ± z)j It is assumed here the color 
charge distribution in x~ , due to the Lorentz contrac- 
tion of the light cone charges, is well approximated by 
a delta function. Also, p a (xT) is the color charge dis- 
tribution of the hard static sources (at large Bjorken x) 
in the transverse plane. These sources are coupled to 
the "soft" (small Bjorken x) partons and different color 
charge configurations are assumed to be distributed with 
the Gaussian weight 



P[p) = exp 



dx d 2 XT 



p a (x ,x T )p a {x ,x T ) 
2A 2 s (x-,x T ) 



(1) 

In the original McLerran-Venugopalan model, the cor- 
relator of color charges in the transverse plane of the 
nucleus satisfied the relation 



(p a (x T )p b (y T ))=A 2 s S^(x T -y T )5 ab 



(2) 



where is a constant. As pointed out by Lam and 
Mahlon |3^], the 5-function form of the uncorrelated 
color charges is inconsistent with color neutrality at 
large transverse separations. In general, we will have 
( P a {x- 7 x T )p b (0)) = A 2 s (x-,x T )6 ab . 

In the limit of an ultrarelativistic source, the classical 
equations of motion for the wee parton fields are 



p a (x ,X T ), 



(3) 



where F^- a = d^A u > a 



d v A^ a + f abc A^ b A v < c is the 
field strength tensor and = + iA^ is the covariant 
derivative associated with the wee parton gauge field A 11 . 



We first find a solution which satisfies A 



c 



A 



c 



in 

the covariant gauge d^A^ = 0. For this solution, the 
covariant gauge condition becomes d+A^, = 0. Then A q 
satisfies the Poisson's equation 



-V T Aj = p(x , xt) ■ 



(4) 



Note that one can alternatively obtain Eq. (Q) in A~ = 
gauge with the static condition d + A^ — 0. Gauge 
transforming our solution from covariant gauge to the 
light-cone gauge A + = 0, one finds the solution [|[ ||, ||, 



A' 



0, A l = iUd l U\ id-U ] = KU\ (5) 



where U is a gauge transformation matrix in the funda- 
mental representation. Here we define A = A q. Also, 
note that in the light-cone gauge, A 1 (i = x, y) is a pure 
gauge in the transverse direction: F xy = 0. The solution 
of the last of the equations in Eq. (S) gives 



U = P exp 



A(z , xt) dz 



(G) 



From Eq. (@), one then has finally the solution 

-1 




2 p(z ,xr)dz~ 



=„p{z ,x T )dz~ 



i\ t 



(7) 



The only non-zero component of the field strength is F +l . 

Thus far we have not imposed any restriction on the 
color charge distribution. If we impose the simple and ob- 
vious constraint that the color charge distribution must 
be zero outside the nucleus, the solution of Poisson's 
equation in Eq. (Q) can still give a non-zero gluon dis- 
tributions outside the nucleus. In two dimensions, the 
fall-off of the gluon field is rather slow as shown in Fig. |^. 
This slow fall-off is a problem for a finite nucleus since the 
gluon field is associated with a non-zero field strength. 
Clearly the simple prescription for color neutrality is not 
sufficiently stringent. 

A more realistic prescription would be to apply the 
color neutrality constraint already at the nucleon level. 
Our numerical procedure to implement the constraint for 
finite nuclei is as follows. We first sample A nucleons on a 
discrete lattice requiring that they satisfy a Woods-Saxon 
nuclear density profile in the transverse plane. Note that 
this procedure generates the same distribution in the con- 
tinuum as 



A 2 s (x T ) = A 2 s0 T A (x T ) 



(8) 



where Ta{xt) = dz n(r) is a thickness function, xt 
is the transverse coordinate vector (the reference frame 
here being the center of the nucleus), n{r) is the Woods- 
Saxon nuclear density profile, and A^ is the color charge 
squared per unit area in the center of each nucleus. The 
only external dimensional variables in the model are A s o 
and the nuclear radius R. 

Next, Gaussian color charge distributions are gener- 
ated on the lattice, where the probability distribution of 
color charge in a nucleon is expressed as 



P[p] = exp - 



Pi 



N 

^2A^~ 

J n -3 



(9) 



where A^ j is the color charge distribution squared, per 
unit area, of a nucleon at a lattice site j and N is the 
number of lattice sites that comprise a nucleon. A^ • is 
obtained from A^g by assuming that the color charges of 
the nucleons add incoherently. There are two versions of 
the subsequent step. In the first (which we term Color 
Neutral I), we subtract from every pj the spatial average 
J2j Pj/N in order to guarantee color neutrality (p) = 
for each nucleon. In the second (termed Color Neutral 
II), the dipole moment d of each nucleon is eliminated 
in a similar manner, by superimposing the result of the 
Color Neutral I procedure with a uniform distribution 
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FIG. 1: Color Charge Correlator in momentum space. Orig- 
inal MV model is shown by circles, while squares correspond 
to the result from color neutrality condition I and triangles 
correspond to color neutrality condition II (see text). The 
color charge correlator is plotted versus pr in units of A s q- 



FIG. 2: Gluon field as a function of radial distance. Original 
MV model is shown by circles, while squares correspond to the 
Color Neutral II prescription (see text). Results from Eq. dlj| ) 
are shown by triangles. The results are for A s qR = 18.5. 



of dipole moments, whose intensity, integrated over the 
nucleon, is — d. This in turn amounts to placing the ap- 
propriate charges along the line boundary of the trans- 
verse cross-section of the nucleon. There is a priori no 
reason to prefer one prescription over the other as long 
as neither leads to significant field strengths outside the 
nucleus. 

In Fig. |l|, we plot the Fourier transform of the charge 
correlator, which in the continuum is defined as 



p(.Pt) 



d 2 x T exp{ip T ■ x T )(p a (x T )p a (0)) , (10) 



for the MV model and for the two variants which impose 
color neutrality on the nucleon level. In the MV model, 
it is evident from Eq. (|J) that p(pr) is a constant every- 
where except at pt = where it is constrained to be zero 
from the global charge constraint. In the Color Neutral 
I (II) variant, we see that p(pr) ~ p\ (~ Pt) f° r sma ll 
momenta pt < A s q and is constant at larger momenta. 
The oscillatory behavior seen for Color Neutral I and II is 
due to the fact that the correlator in coordinate space is 
not strictly a delta-function. In the coordinate space the 
charge correlator for the two models, Color Neutral I and 
II, falls off rapidly, as ~ l/x T and ~ l/x T respectively, 
at larger distances. 

It is an interesting coincidence that the behavior 
of p{pt) in our model is similar to the behavior ex- 
pected from the renormalization group (RG) evolution 
of color charges in the McLerran-Venugopalan model. In 
Ref . p2l , it is shown that the screening of color charges 
due to the RG evolution gives a behavior p(pr) ~ p\ 
for pt < A s o (and /5(;>r)=constant for pt > A s o). A 
recent analysis by Mueller gives the same result with a 



specific prediction for the prefactor [Q . These results lie 
between the range spanned by our "Color Neutral I" and 
"Color Neutral II" . Since the only information on the RG 
in our formalism for nuclear collisions comes from the ini- 
tial conditions, our results may be quantitatively similar 
to RG evolved predictions for nuclear collisions. 

We now discuss how the gauge fields are determined 
from the color charge distributions. The pure gauges at 
the site j are defined with Hermitian traceless matrices 
Aj on the lattice as 



Un.j = exp[iA 



j+n\ exp[ 



-iAi 



(11) 



Aj can be obtained by solving the lattice Poisson equa- 
tion 

-AiAj = -^2 (A J+ „ + Aj_„ - 2Aj) = Pj . (12) 

n— x ,y 

We solve Eq. (|l2| ) numerically, using the fast Fourier 
transform (FFT) method which is faster than the ap- 
proach employed in Ref. fl5j| . Imposing color neutrality 
on the nucleon level (for both Color Neutral I and Color 
Neutral II) significantly suppresses the gluon field outside 
the color charge density as shown in Fig ||. 

One can also simulate our color neutrality prescription 
by introducing a screening mass m c to the Poisson equa- 
tion Eq. (H) 



(-A L +ml)A j = p j , 



(13) 



as a simple method to suppress the gluon field. The 
range of the screening mass must be m c ~ Aqcd- Al- 
though this procedure does not originate in any color- 
neutrality condition, it is clear that the gluon distribu- 
tion at high momentum range should not be modified, 
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FIG. 3: Comparison of the unintegrated gluon distribution 
before the collision obtained by the simple MV model (tri- 
angles), color neutral I (circles) and color neutral II (full 
squares). The SU(2) result scaled by 16/9 is shown in open 
squares. The line corresponds to a 1/pj. distribution. Results 
are for A a oR = 83.7. 



since |At(pt)| 2 ~ l/(PT + m c)- Fig || demonstrates that 
the gluon field distribution with m c = 1.0 fm -1 is nearly 
identical to the color neutral distribution. 

In Fig. ||, we plot numerical results for the SU(3) un- 
integrated gluon transverse momentum distributions be- 
fore the collision. The gluon distributions from models 
which impose strict color neutrality at the nucleon level 
are significantly suppressed at momenta pt < A s o com- 
pared to that of the original MV model which had only a 
global color neutrality constraint for the nucleus. In the 
large momentum region, where we expect the color neu- 
trality constraints to be immaterial, both sorts of models 
show the same 1/px dependence. Comparing the SU(2) 
result for the Color Neutral II model to the SU(3) result, 
we observe that the ratio of the two results is the ratio of 
the respective Casimirs. The results shown in the figure 
are all for one value of the dimensionless coupling A s oi?. 
The qualitative results arc the same for other values of 
this parameter. 



B. Nuclear collisions 

We can now apply our model for parton distributions 
of finite ultrarelativistic nuclei to discuss classical gluon 
production in nuclear collisions. This problem was first 
formulated by Kovner, McLerran and Weigert |TlJ and 
perturbative solutions were discussed further by several 
authors [jl2], |l3|, The numerical formulation of the 
classical problem has been discussed extensively by us 
in previous papers [H| [l(| [l?], Q. In these papers, we 
considered, for simplicity, only the case of cylindrical nu- 
clei with uniform nuclear matter distributions. Periodic 



boundary conditions were imposed. We will discuss here 
the formulation of the problem for finite nuclei and open 
boundary conditions. In the previous sub-section, we ob- 
tained the color charge and gluon distributions in coor- 
dinate space and momentum-space after applying color 
neutrality prescriptions at the nucleon level. These distri- 
butions then give us the necessary ingredients to discuss 
the space-time evolution of classical fields after the col- 
lision. For completeness (and continuity) we will repeat 
some of the arguments in the papers cited here. 

We fix the gauge condition to be the radiation gauge 
A T = or equivalcntly, 



c + A~ 



'A+ = 0. 



(14) 



We require that the light cone color charges on each light 
cone are delta function sources: J ±,a = 5(x : ^)p ± ' a (xT)- 
the solutions to the Yang- Mills equations are then ex- 
plicitly boost invariant. We are interested in the physics 
pertaining to the mid-rapidity region and will consider 
the space-time rapidity r\ to be zero. Our ansatz for the 
gauge field as a function of proper time r = \f2x + x~ is 



.4' 



A ± = 



4(t, x t )6(x-)9(x+) 
a\(T,x T )6(x-)6(-x + ) 
ai 2 {T,x T )0(-x-)9(x + ), 
±x ± a{r, x T )0{x-)0{x + ). 



(15a) 
(15b) 



The boundary conditions are determined by matching the 
solutions in the space-like and time-like regions. Requir- 
ing that the gauge fields must both be regular at r = 0, 



Oand D^+F^+ = J+ for x~ 
boundary conditions at r = 0: 



gives the 



4(0, xt) 
a(0,x T ) 



a\(0,x T ) + ai(0,XT), 



(16a) 



[a\(0,x T ),ai(0,x T )]. (16b) 



Note that these conditions, first formulated for infinitely 
large nuclei, are the same for finite nuclei. There is no 
ambiguity in determining the reaction zone since, for ex- 



ample, when a\ — 0, we have a — 0, 



a 2 , a pure 



gauge-hence there is no classical gluon production out- 
side the reaction zone. 

We shall now explain how one generates the reaction 
zone in the model. The initial conditions (at r = 0) can 
be generated in two possible ways. 

(a) For a given impact parameter 6, one determines the 
number of participant nucleons using the Glauber 
model as for example in Ref . |^] . The corresponding 
color charge density (and the gluon distribution) 
for each nucleus in the reaction zone can then be 
determined using the sampling procedure described 
previously. The matching condition in Eq. ([l6]) is 
then used to determine the initial conditions for the 
evolution. 

(b) The color charge distribution is determined sep- 
arately for each nucleus by projecting the color 
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FIG. 4: The two components of the transverse pressure (T xx 
and T yy ) and the energy density e plotted as a function of r in 
dimensionless units. The results are for a impact parameter 
b — R and dimensionless coupling A s oR = 18.5. Also shown 
is the sum of the two transverse pressures. 



charge distributions of all the color neutral nucle- 
ons in the transverse plane. The gluon distribu- 
tion for each nucleus is then determined accord- 
ingly, and the reaction zone is automatically given 
by the matching condition in Eq. ( fl6| ) for a speci- 
fied impact parameter b. 

The reaction zone can be determined naturally from 
the procedure (b) as long as we have the correct color 
charge distribution for a finite nucleus, since the region of 
A\ = 0, A 2 remains pure gauge. Nevertheless, we would 
like to use (a) for purely technical reasons, namely, since 
it is more efficient for numerical computations. We have 
checked that the two procedures give nearly identical re- 
sults. 

Explicit expressions for the discretized lattice Hamil- 
tonian, for other components of the stress-energy tensor, 
and the classical equations of motion can be found in the 
appendix. In Fig. |], we plot the temporal behavior of 
the two transverse pressures and the energy density. In- 



terestingly, at late times we notice that e = T xx + T % 



as one would anticipate for a free streaming gas of glu- 
ons in the transverse plane. We will return to this result 
later. The difference between the two transverse pres- 



Tr 



Tyy is a measure of the azimuthal 



anisotropy generated by classical fields pOl 



III. IMPACT PARAMETER DEPENDENCE OF 
MULTIPLICITY DISTRIBUTIONS 

In previous papers JlTJ, [Tt], [l^], we compared our re- 
sults for cylindrical nuclei with periodic boundary con- 
ditions and uniform nuclear matter distributions to the 
RHIC data. In this section, we will apply the more re- 



TABLE I: A s0 = 1.41 GeV. Results are for a 256 x 256 lattice- 
the nuclear radius is 64 lattice units. All dimensionful scales 
are in GeV units unless otherwise stated. 



b(fm) 


N part 


g 2 N g 


g 2 E g 


A(6) 


Qs(b) 


Mb) 


0.000 


377.89 


1628.68 


2725.40 


1.1777 


1.0836 


0.3695 


3.150 


321.35 


1309.83 


2170.08 


1.1545 


1.0563 


0.3484 


4.725 


263.33 


1035.84 


1663.88 


1.1234 


1.0197 


0.3359 


6.300 


199.11 


760.95 


1182.02 


1.0741 


0.9623 


0.3249 


7.875 


136.47 


515.61 


751.902 


0.9993 


0.8758 


0.3152 


9.450 


81.21 


295.76 


384.004 


0.8876 


0.7487 


0.3004 



TABLE II: A s0 = 2.32 GeV. Results are for a 512 x 512 lattice- 
the nuclear radius is 128 in lattice units. All dimensionful 
scales are in GeV units unless otherwise stated. 



b(fm) N part g 2 N g g 2 E g 



A(b) Q„(b) f N (b) 



0.000 
3.150 
6.300 
7.875 
8.367 
9.450 



377.89 3768.00 
321.35 3061.59 



199.11 
136.47 
118.17 
81.21 



1808.89 
1215.17 
1042.54 
699.95 



9198.736 1.9517 2.0355 0.2867 

7492.084 1.9132 1.9866 0.2757 

4183.888 1.7800 1.8185 0.2610 

2636.300 1.6560 1.6636 0.2522 

2243.692 1.6060 1.6017 0.2514 

1411.900 1.4708 1.4356 0.2411 



alistic framework developed in the previous section to 
compute the gluon multiplicity distributions as a func- 
tion of centrality. If one assumes parton-hadron duality, 
these distributions can be compared, for instance, to the 
centrality dependence of the charged hadron spectra in 
heavy ion collisions at RHIC ||. 

At late times r ^> l/A s o, the dynamics of the classical 
Yang-Mills fields produced in the collision can be lin- 
earized and approximated by that of a system of weakly 
coupled harmonic oscillators. In previous papers |l7|, p"8] , 
we discussed two methods to compute the gluon number. 
Both of these are applicable when the harmonic oscilla- 
tor approximation is a good one. One method is a gauge 
invariant relaxation method which can only be used to 
measure the particle number and not its distribution with 
momenta. The other method is to fix the Coulomb gauge 
and compute the field amplitudes squared in momentum 
space. We showed that both methods give results for the 
gluon number which agree very well with each other. 

We will therefore restrict ourselves in this paper to 
using one method alone-that of fixing the transverse 
Coulomb gauge V_l ■ A = [0 and computing the field 
amplitudes. We have 



N(k) = u,(fc)Mfc)| 2 ) = vfOTPH, (17) 

where <fi(k) and ir(k) correspond to the potential term 
and the kinetic term in the Hamiltonian respectively. 
The average () represents an average over the initial con- 
ditions. 

In Refs. 0, [l8|], non-perturbative relations were de- 
rived for the energy and number of produced gluons (at 
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central rapidities) as a function of A^. For realistic nuclei, 
these non-perturbative relations are less simple. One can 
parametrize our results for the gluon number with the 
more general relation 



dNg 

dr) 



f N (b) / d 2 X T 



9 2 ' 



(18) 



where A s (b, xt) is the local saturation scale defined to 
be A 2 (b, xt) — C ■ p(b,XT)/2, where p is the participant 
density at a particular position in the transverse plane, 
and C is the color charge squared per nucleon. When 
A s (b, XT)=constant, as for cylindrical uniform nuclei, one 
recovers the form of the expressions in Refs. |17, Il8). 



The color charge squared in the center of the nucleus is 
A 2 s0 = C ■ p(0,0)/2 : so A 2 s (b,x T ) = A 2 s0 p(b,x T )/p(0,Q). 
One can then re-write the previous equation as 



dN g _ f N (b) A% N 



d-q 

p(0,0) 



Po 



part 



(b), 



(19) 



4.321 fni and N, 



part 



where po = 
J d 2 x T p(b, x T ) 

In Tables [| and y, we show the calculated SU(3) re- 
sults for two values of the saturation scale in the center 
of the nucleus: A s0 = 1.41 and A s0 = 2.32 GeV respec- 
tively. In the tables, b is an impact parameter (in units of 
fm) and N part is a number of participants at that impact 
parameter. The latter is calculated using a Woods-Saxon 
nuclear density profile. We list in the tables our results, 
as a function of impact parameter, for g 2 N g ; the number 
of produced gluons and g 2 E g ; the transverse energy of 
produced gluons in GeV multiplied by the value of the 
strong coupling constant squared g 2 , evaluated (to one 
loop order) at the average value of the saturation scale 
(denoted in the tables as Q s (b)) for that impact param- 
eter. 

In Table |, we see that the gluon multiplicity, for the 
most central collisions, is about half the final pion mul- 
tiplicity at RHIC while the transverse energy of gluons 
is close to the experimentally observed transverse energy 
of about 600 GeV. In Table ||, the gluon multiplicity is 
close to the final pion multiplicity at RHIC while now 
the transverse energy is nearly four times larger than the 
experimentally observed transverse energy. Neither the 
initial transverse energy nor the initial gluon multiplic- 
ity are conserved quantities, so the results for different 
values of the saturation scale lend themselves to very 
different dynamical scenarios. For instance, the authors 
of Ref. [g0| have suggested that inelastic gg — > ggg pro- 
cesses may be important in the early stages of nuclear col- 
lisions thereby increasing the gluon multiplicity prior to 
thermalization/hadronization. Alternatively, if the gluon 
number close to the experimental number (as in Table ||) , 
strong hydrodynamic expansion at early times might be 
essential to lower the initial transverse energy to the ex- 
perimentally desired value. Additional observables (such 
as v 2 ) may be important in distinguishing between the 
different scenarios. 
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FIG. 5: Comparison of the centrality dependence of the gluon 
distribution from SU(3) lattice results to data from experi- 
ments 0, [35) . The strong coupling constant is fixed to the 
value g = 4. The lattice results for A s (0) = 1.18 GeV and 
A s (0) = 1.95 GeV are multiplied by a factor 2.4 and 1.1, 
respectively. 



In Fig H, we plot the centrality dependence of gluons 
calculated using Eq. (17) together with the experimen- 

11 and PHENIX f 



tal data from PHOBOS |4J and PHENIX |35|. We as- 
sume here that the charged particle multiplicity is two 
thirds of the gluon number. The classical computation 
is performed for fixed a s ; the centrality dependence, as 
seen from Eq. @, comes from the dependence of /at on 
the impact parameter. In Ref. [fl7|| , it was shown that 
/jv = /jv(A s i2) increases slowly with A s i?-hence one ex- 
pects it to vary with impact parameter. We see that the 
results agree reasonably well with the data. 

The centrality dependence of the transverse energy is 
studied in Fig. |^. As in the case of the multiplicity, even 
though the absolute normalization is strongly dependent 
on one's choice of A s , the centrality dependence is very 
similar for the two A s 's and shows reasonable agreement 
with the data. 

Let us now compare our results with those derived in 
Ref. [2(j] and discussed further in Ref. [0. In these 
works, one obtains in terms of Q s , the average satura- 
tion scale, the result 



dr/ 



N 2 -l f, 2 Q 2 s (b,x T ) 
4n 2 N c J XT a s 

/ ^9/j\ \ ^ Dart 



CN xG(x,Q 2 s (b))^. 



(20) 



In the leading logarithmic approximation, if cat is con- 
stant, one obtains a logarithmic dependence on the cen- 
trality entirely from xG(x,Q 2 (b)). One could thus at- 
tribute the logarithmic behavior at the classical level to 
fixed a s and leading logarithmic behavior of the gluon 
distribution function or equivalently, at higher order, to 
the one loop running of a s . In the physically interesting 
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FIG. 6: Comparison of the centrality dependence of the gluon 
transverse energy distribution from SU(3) lattice results to 
data from experiments [^6). The strong coupling constant is 
fixed to the value g 2 = 4. The lattice results for A 3 (0) = 1.95 



GeV are scaled by j 

by I 



while those for A s (0) = 1.18 are scaled 



regime, it is difficult to distinguish between the two. A 
reasonable agreement with the data is also seen in this 
formulation of the problem. 

The relation between the two formulations is as follows. 
As discussed previously, what we simulate numerically 
is the color charge squared per unit area, A 2 fl7fl . The 
saturation scale, on the other hand, is a scale determining 
the behaviour of the gluon number distribution ]29[ ] . The 
two scales can be related by computing, for instance, the 
gluon number distribution in a nucleus in the McLerran- 
Venugopalan approach (where A s naturally appears) and 
in Mueller's dipole approach. On comparing the two one 
obtains, to leading logarithmic accuracy, the following 
relation between A s and Q s : 



Qt 



Ait 



log 



A? 



A 2 

QCD 



(21) 



with Ar 



^QCD — 0.2 GeV. The relation between c/v and Jn 
is c N = An 2 f N /{N 2 c - l)/ln(Q 2 /A| CD ). Therefore, if 
cj\r is to be a constant, fpf increases logarithmically with 
Q s . A weak rise in /jy is seen in our simulations. If the 
infrared scale, as argued in Ref. j}7|, is a number of order 
0(1/Q S ), we would have A s ss Q s . Note that in practice 
Q 2 s {b) and A 2 (&) (related by Eq. ( |2l| ) are nearly equal to 
each other for all impact parameters. 

In Ref. [pp[ , it was argued that subsequent inelastic 
scattering of gluons after the classical approximation is 
applicable but before thermalization would increase the 

2/5 

multiplicity by a factor l/a s ■ This effect was taken 
into account in Ref. |^lj and it was shown that good fits 
to the data could be obtained for this functional form as 
well. The effects of parton multiplication a la Ref. plf] 
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FIG. 7: Top: Gluon transverse momentum distributions as a 
function of p T /A s (b) for b = 0.0, 6.4, and 9.45 fm. Bottom: 
Same as top figure but with re-scaled normalizations. 



are beyond the scope of this work and we will not discuss 
it further here. 

We now turn to the computation of the initial gluon 
transverse momentum distributions. In the top half of 
Fig. f^, we show the initial gluon transverse momen- 
tum distributions for different impact parameters in di- 
mensionless units of pr/A s (b). In the bottom half of 
the figure we plot the same quantity by re-scaling the 
overall normalization of the distributions. Strikingly, as 
a consequence of geometrical scaling, the distributions 
now coincide-the initial gluon distributions are a uni- 
versal function of pr/A s (b). In Ref. ^7j, it was argued 
that the RHIC data display such a geometrical scaling. 
More detailed comparisons with the recent RHIC data at 
y/SNN = 200 GeV will help further clarify the observa- 
tion of geometrical scaling. Note however that the typical 
(px) of the initial gluon distribution is of order A s , while 
the observed (pt) of charged hadrons is smaller, of the 
order 3T C . For the A s in Table |, this difference is a factor 
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FIG. 8: Top: Gluon transverse momentum distributions at 
r\ — for impact parameters 6 = 0, 6.3, 9.45 fm in the case of 
A s0 = 1.15 GeV. Mini-jet partons from PYTHIA, multiplied 
by number of collisions, with cut-off momentum pr = 2.2 
GeV/c, are also plotted in histograms. Bottom: The same 
for impact parameters b = fm for A s o = 1-15 GeV and 
A 3 o = 2.23 GeV. Different cut-off momenta of pr = 2.2 and 
Pt = 1-8 GeV/c for the PYTHIA mini-jets are also plotted 
by histograms. 



of two; for those in Table ffl it is significantly larger. 

It is instructive to compare the transverse momentum 
distributions computed in the classical approach with the 
pQCD computations of the gluon distribution. The latter 
is obtained as follows. In the parton model, when the 
parton density is small, the cross section for hard parton 
scattering can be calculated as 



da 



je.t 



dp T dy\dy 2 



u dt 



(22) 

where y\ and y 2 are the rapidities of the scattered par- 
tons and X\ and x 2 are the fractions of momentum of 
the initial partons. The structure functions, f a (x,Q 2 ) 



are taken to be the CTEQ5 leading order parametriza- 
tions The parton distributions are evaluated at 

the scale Q 2 = p T . Here da a b/dt is the leading order 
QCD partonic cross section. In addition, one includes 
initial and final state radiation processes in order to take 
into account the enhancement of higher-order contribu- 
tions associated with multiple small-angle parton emis- 
sion. This is for instance what is done in the event gen- 
erator PYTHIA and we use for our purposes PYTHIA 
version 6.2 |3!j for the incident cm. energy of ^/s = 130 
GeV. The K factor of K = 2 is chosen to fit the UA1 
data of pp at y/s = 200 GeV @. The nuclear gluon 
distributions are obtained by multiplying the PYTHIA 
result by the number of hard binary collisions. Nuclear 
shadowing effects are neglected. 

In the top half of Fig. |^, we study the impact parameter 
dependence of the classical gluon transverse momentum 
distributions (with A s o = 1.14 GeV) together with the 
parton distributions from PYTHIA calculations with an 
infrared cut-off parameter of pr = 2.2 GeV/c. We find 
that there is a overlap region between 1-3 GeV/c between 
the two computations. For higher transverse momenta, 
the classical computation overshoots the PYTHIA re- 
sult. One reason this is the case is that the realistic 
parton distributions have not been included in the clas- 
sical distribution-this makes a difference at large trans- 
verse momenta g§ ||. Naively, one would expect the 
pQCD result to diverge rapidly at small px- The reason 
it doesnt is that initial and final state radiation processes 
change the parton distributions from Eq. ( |22] ) below 3-4 
GeV/c. In the lower half of Fig. ||, we plot the gluon dis- 
tributions with the parameters A s o = 2.32 and 1.14 GeV 
along with the PYTHIA gluon distributions for different 
values of the infrared cut-off. The latter is not an infrared 
safe quantity so it is not too surprising that the shape and 
magnitude of the distribution does depend on the cut-off. 
So if one were to match the classical and pQCD distribu- 
tions, one has to adjust the px cut-off for different values 
of A s o- The differences in the gluon number between the 
classical and classical+PYTHIA LO pQCD result are not 
very significant for the particle number but may be so 
for the transverse energy since the latter is more ultravi- 
olet sensitive. Therefore the large transverse energy per 
particle obtained in the classical approach may be sig- 
nificantly overestimated. A more realistic calculation of 
the distribution of hard modes will find less of an excess 
(relative to experiment) in the transverse energy. 

We discussed in this section our results from numerical 
simulations of very high energy nuclear collisions of nuclei 
at various impact parameters. While the absolute num- 
bers from our results depend very sensitively on the scale 
A s , the dependence of the gluon number and transverse 
energy on centrality and the phenomenon of geometrical 
scaling of the momentum distributions are independent 
of A s and appear universal. In addition, they appear 
to agree reasonably well with trends in the RHIC data. 
Similar observations were made in Refs. |)[ |26| and in 
Ref. [p7| for global features of the data. 
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IV. PROTON-NUCLEUS COLLISIONS 

Collisions of protons/deuterium with large nuclei at 
RHIC energies are essential to distinguish dynamical ef- 
fects due to correlations in the nuclear wavefunction from 
dynamical effects resulting from the collective expan- 
sion of hot and dense matter in nucleus- nucleus colli- 
sions. The first computation, in the classical framework, 
of gluon production was performed by Kovchegov and 
Mueller p9| . Here we will follow the subsequent work of 
Dumitru and McLerran fl30| which is closer in its formu- 
lation to the discussion here. In Rcf. |3(|, the authors 
introduce two saturation scales which we will denote as 
A s i (for the proton) and A S 2 for the nucleus. One ex- 
pects that since the parton density in a nucleus is larger 
than in a proton, <C A^ 2 . In the kinematic region 
Pt > A s i, analytical computations of gluon production 
are feasible. For A s i < A S 2 < pt, the situation is anal- 
ogous to nucleus-nucleus collisions and the distributions 
are proportional to A^A^/p^. For A^i < px < A s2 , 
the distribution changes qualitatively and is proportional 
to A. 2 s1 /pt- In this kinematic region, one obtains ]4l| ] 
dN/d 2 bd 2 k t oc ln(kf/A 2 sl )/k 2 . 

In our framework, pA collisions can be simulated 
straightforwardly since the only change from nucleus- 
nucleus collisions is to consider asymmetric values of A s . 
The results of our simulations are shown in Fig. 0. We 
choose a fixed value for A s o2 in the center of the nucleus 
and plot the momentum distributions for different A s oi's. 
At largely, as anticipated, the distribution goes as l/p T - 
At smaller px of pt ~ 2 A s o2 one sees a qualitative change 
in the distributions which is more consistent with al/p 2 -. 
behavior of the distributions. Thus it appears that our 
distributions reproduce the shapes predicted by Dumitru 
and McLerran in the appropriate kinematical regions. 

A comparison of absolute normalizations is more 
tricky. If one computes analytically j|l| the ratio of the 
multiplicities at the scale k t ~ A s o2 for the two values of 
A s oi shown in Fig. || the ratio is ~ 23, which is the scal- 
ing factor used to scale the distribution with the lower 
A s oi to that with the larger value. However, we should 
caution that the absolute normalization of these results 
is uncertain because the behavior in the kinematic region 
Pt < Asoi cannot be computed analytically jl^] and can 
only be computed numerically as discussed here. 

More detailed studies of pA collisions are underway 
and will be reported separately E3| . 

V. SUMMARY 

In this paper, we extended previous results for colli- 
sions of uniform, cylindrical nuclei with periodic bound- 
ary conditions to finite nuclei with realistic nuclear mat- 
ter distributions and open boundary conditions. To 
achieve this, we imposed color neutrality constraints on 
nucleon color charge distributions. With these improve- 
ments, we were able to study nuclear collisions at differ- 




FIG. 9: Transverse momentum distribution of produced glu- 
ons in a pA collision. Here A^ 01 refers to the color charge 
squared per unit area in the center of the proton while A502 
is the color charge squared per unit area in the center of a 
nucleus. 



ent impact parameter and proton-nucleus collisions. We 
obtained results for the initial gluon number and trans- 
verse energy and compared these to the data. While the 
absolute normalizations of these quantities depend sen- 
sitively on the saturation scale A s , their centrality de- 
pendence appears universal and shows reasonable agree- 
ment with the experimental data. Similarly, the trans- 
verse momentum distributions display geometrical scal- 
ing. (The latter has been argued to be a feature of the 
RHIC data ^7|.) Our results however also suggest that 
dynamical effects at later times (possibly beyond the va- 
lidity of the classical approach) might also be important 
in describing the RHIC data in detail. This is suggested 
for instance by the necessity of reproducing both the 
transverse energy and the multiplicity but more emphat- 
ically by the failure of the classical approach to explain 
the elliptic flow data fiofl . It has been argued that the 
elliptic flow data may be explained by non-flow correla- 
tions |E3 but more detailed theoretical and experimental 
analyses are required to clarify the issue. 

We also briefly discussed proton-nucleus collisions in 
the classical framework. We showed that our numerical 
simulations agree qualitatively with previous analytical 
work |2^, ^0). This agreement suggests that results of 
our numerical simulations for nucleus-nucleus collisions 
(where analytical computations are more difficult) are 
robust. 

The classical distributions discussed here can be used 
as the initial conditions for the subsequent time evolution 
of produced gluons as might be described, for instance, 
in Boltzmann type transport calculations [Q, [2(| |||. HJ] 
at late times where held strengths are weak. With these 
initial conditions, it will be interesting to study the con- 
sequences of parton evolution on flow or parton energy 
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loss. 

The classical treatment of nuclear collisions can also be 
improved significantly. One such improvement is to give 
up the assumption of strict boost invariance and to study 
the effects of doing so on our distributions. Another is to 
go beyond the simple assumption of Gaussian-distributed 
color charges and to incorporate the effects of renormal- 
ization group evolution of the color charges at high ener- 
gies. All of these studies, tempered by the large amount 
of available experimental data, will hopefully provide in- 
sight into the very earliest moments of a high energy 
nuclear collision and of the subsequent dynamics of the 
produced quark gluon matter. 
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APPENDIX: SPACE-TIME EVOLUTION OF THE 
STRESS-ENERGY TENSOR IN A NUCLEAR 
COLLISION 

For completeness, we provide in this appendix explicit 
expressions for the lattice Hamiltonian, the lattice equa- 
tions of motion and discretized expressions for other com- 
ponents of the Stress-Energy tensor for the problem of 
interest. 



density in A T = gauge is 

C = Tr ^-{d TAl f + ±-{d T A v f - I** - i-f*) , 

(A.2) 

where i,j runs over transverse coordinate x and y. 

Now let us assume 77 independence of the fields. As 
discussed previously, the Yang-Mills equations have this 
property if the sources are strictly <5-function sources on 
the light cone. We have 

MT,V,r) = Ai(r,Oi A,(r, V,r) = *(r,f) , (A.3) 

resulting in F vi — —Di<&, where Z>; = di — ig[A i: ■ ■ ■] 
is the covariant derivative in the adjoint representation. 
Defining the conjugate momenta 

Ei = Td T A h Pv = -d T A„ , (A.4) 

T 

one finds that the boost invariant Yang-Mills Hamilto- 
nian is the QCD Hamiltonian in 2+1 dimensions coupled 
to an adjoint scalar p"5[ : 

H=jd^{^Ef + ^ + T -Fi j + ^rY 

(A.5) 

The equations of motion corresponding to this Hamilto- 
nian are 

E 

M = — , $ = rp V) (i = x,y), 

T 

E x = -TD v F xy + -[$,D x $}, 

T 

E y = rD x F xy + — [$, D y <f>], 
Pv = l(D x $ + Dl$), 

(A.6) 

and the Gauss law constraint is 

D-E + ig\pr,,$] =0. (A.7) 



1. Hamiltonian approach in the transverse plane 

In this sub-section, we write down the classical Hamil- 
tonian's equations of motion in the r, 77, r coordinates 
where r = ^/2x + x~ is the proper time, f] = \ log(x + /a;~) 
is the space-time rapidity and r = (x, y) are two trans- 
verse Euclidean coordinates. Inverse transformations are 
t = r cosh 7/ and z = rsinh?7 respectively. The QCD 
action for gauge field in the r, 77, r coordinates reads 



Tdrjdrd 2 



1 



-Tr(g^g^F af3 F^) 



(A.l) 



where F^ u = d^A v — d v A^ — ig[A lll A v \ and the metric is 
diagonal with g TT = —g xx = —g yv — 1 and g m = —1/t 2 . 
A^ = A^ t a , and t a represent a gauge group matrices with 
the normalization of Tr(t a t b ) — 2&y. The Lagrangian 



2. Discretized Hamiltonion formalism in 2+1 
dimensions 

In order to realize numerically the solutions to the 
equations of motion in the previous section, while main- 
taining the gauge symmetry, we introduce the link vari- 
ables at the site i 

Uj,i = exp [igaAj(i)] , (j = x, y), (A.8) 

where, a is a lattice spacing. Defining the plaquette 

Un^U U] U mj+l ul J+m Ul^ (A.9) 

the magnetic energy in the transverse plane is expressed 
for a — ► as 



2aV ^ 



(N c - ReTrU c 



(A.10) 
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where the relation Uu — exp[ia 2 gF^,y] was used. The 
Hamiltonian on the lattice is 



H L = ±- £ Tr^ 



-^(iV c -ReTr?7 n ) 

□ 

^Tr($ 3 -f/ Aj <D 3+A C/] J f 



E T ^ 2 > 



(A.ll) 



where the following dimensionless variables are used |4q, 



H - A. ^ ^> ( A - 12 ) 
$ — ► — , r — > flT. 

Hamilton's equations of motion on the lattice are then 

i 

U r ,i = -~U rti E ri i, ®i = Tpi, (r = x,y), 



x.i+y y,i 



Ul m iUl : i-yUx,i-yU yj i+x-y~^ 



"I" U^ i+ Q_ S .Uyj_.Ux,i-£Uy > i\ 

Pi = -[ E ( u A4®*+i u b + u ki-n®i-* u ^-*) 



4$, 



(A.13) 



Once one has explicit expressions for the fields and 
their canonically conjugate momenta, other components 
of the Stress-Energy tensor can be computed as well. 
In light-cone coordinates (r,r],r), the symmetric energy- 
momentum tensor can be defined as 



= Tr(F M " F v a ) — g^ v C . 



(A.14) 



The diagonal spatial components T xx and T yy (the two 
components of the transverse pressure) are explicitly 

T xx = ^Ti(^(-El+E 2 y + (D X $) 2 - (Z^) 2 ) 

+ Tpt + lFijF^, (A.15) 

T yy = lTr(~(+E 2 x -E*-(D x $) 2 + (D y $) 2 ) 

+ rfi + ^FijFij). (A.16) 



The T TX and T Ty parts can be calculated as 



T rx = T^EyF^+PniDz*)), 



(A.17) 



T Ty is obtained by replacing x <-> y in this expression. 

The lattice equivalents of these continuum expressions 
are as follows. The transverse pressure in the y direction 
is 



uu 



I E Kir 



j=i^—x,i—y 



+ ^(iV c -ReTrC/ n ) 

4D 

j—i^i—x 



3=1,1— y 



(A.18) 



where 4D means a sum over the four plaquettes which 
contain the link i with the same orientations. T xx can 



be obtained by replacing x <-> y in Eq. (A.l 



One of the ways to get the lattice definition of Poynting 
vector may be a use of the conservation law d^T^ = 0: 



d T T TT = d x T TX + d v T Ty , 



(A.19) 



with the T TT which has to be defined in a local gauge 
invariant form at a site i as 



rpTT 



8r 



l - ^Ej + T -{N c -KeTrUu) 

sT e=(j,n) 



4r 

r 
16 



3,n 



E Trp j 2 ' 



(A.20) 



where = i + (i + x) + (i + y). After quite some algebra, 
we obtain 



1 



rE, 



+ z[pi(Uli- a! *i-x-Ux,i$i) 
+ Pi+v( U l,i+y- x ®i+V-x 

-U Xi i +y $ i+y+x ) . (A.21) 



U x ,i+yU y!i+x U Xti Uy t ; 
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